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1 Introduction 

First idealized models describing the normal 
activation sequence in the human heart were 
reported over two decades ago [1,2]. Mainly due to 
computational limitations, the models did not 
include myocardial anisotropy nor physiological 
propagation. On the basis of anisotropic bidomain 
theory [3] and cellular automata theory [4], 
development of more realistic whole-heart models 
have become feasible [5]. 

A ventricular model that produces a correct normal 
activation sequence is a prerequisite for simulating 
pathological conditions, such as ischemia, infarction 
or ventricular arrhythmias. A comprehensive model 
should feature an anatomically accurate geometry, 
intramural fibrous structure, and a conduction 
system. Implementation of the ventricular 
conduction system is a challenge, since it should be 
flexible enough to allow “rewiring” of the 
conduction system for each individual case where 
endocardial mapping data are available. Advances in 
computer technology have only recently made this 
approach possible. 

Detailed studies on the anatomy of the human 
conduction system are scarce, and often report a 
high degree of interindividual variability. Therefore, 
we followed the criterion that the conduction system 
should rather accurately reproduce correct the 
pattern of ventricular activation sequence than to 
follow any predetermined anatomical design. The 
conduction system model was validated, first by 
comparing the simulated activation isochrones to 
those reported in the literature, and second by 
comparing the resulting body surface potential maps 
and magnetocardiographic distributions to our 
recordings obtained in normal subjects. 

2 Methods 

2.1 Anisotropic propagation model 

Our hybrid model of the anisotropic myocardium 
describes the subthreshold behaviour of the elements 
according to the bidomain theory [3], while in the 


suprathreshold region the elements behave as 
cellular automata [4]. Each element is assigned a 
specific type and a vector of local fiber direction [5, 
6]. Each element can be in one of four macrostates, 
corresponding to different phases of an action 
potential: resting, excitatory, absolute refractory, or 
relative refractory. During simulation, elements 
undergo a series of macrostate transitions, which 
have been described in detail in Ref. [7]. The 
subthreshold transmembrane potential, v m is 
calculated from the generalized cable equation, 
derived from the anisotropic bidomain theory [3]: 

c m 5v m / 3t = k/ (k+1) V»Dj Vv m - i ion + i app (1) 

Here we assume that the assumption of equal 
anisotropy ratios is valid [5, 7]. In Eq. 1, c m is the 
specific membrane capacitance, Di is the 
intracellular conductivity tensor, and i ion and i app are, 
respectively, ionic and applied currents. Above a 
threshold value, the cell enters the absolute 
refractory state and v m is estimated from a 
predefined logarithmic function [7]. Physiological 
parameters for the model are adopted from literature 
[8,9]. 

2.2 Extracardiac field calculation 

In the forward computations, the ventricles are 
assumed to be embedded in a torso-shaped 
piecewise homogeneous volume conductor, 
including the lungs and the intracardiac blood 
masses. The infinite medium potential, (J) 0 , in the 
extracardiac region is determined from the 
discretized equation 

-47TG(t)o = Gij Vv m • r / r 3 dV 

+ c 2 J a a T Vv m • r / r 3 dV (2) 

where the integrals are evaluated over the 
ventricular volume, v m is as in Eq. (1), a is the local 
direction of the fiber axis, G is the extracellular 
conductivity (g = 2.0 mS/cm), Gi is the transverse 
conductivity (Gi = 0.5 mS/cm) and g 2 is the 
difference between the axial and transverse 



conductivity (o 2 = 1.5 mS/cm) [7]. The first term of 
the potential in Eq. (2) represents the contribution of 
the isotropic component and the second term 
accounts for the anisotropic properties of the cardiac 
tissue. 

The extracardiac magnetic field, H 0 , is evaluated 
from a corresponding equation, 

-47TgH 0 = GiJ Vv m X r / r 3 dV 

+ g 2 J a a T Vv m xr / r 3 dV (3) 

To compute the body surface potential and magnetic 
field in the inhomogeneous torso, a "fast forward 
solution” is used [7]: 

^ = (I-n)- 1 ^> 0 , H = H 0 + A(I-n)- 1 ^> 0 (4) 

In Eq. (4), Q is a geometry matrix consisting of the 
solid angles subtended by a surface element at each 
node, and A is another (vector-value) geometry 
matrix taking into account the influence of the 
conductivity change surfaces on the magnetic field 
(for details, see Ref. [10]). 

2.3 Fiber arrangement 

The model consists of 2.000.000 excitable elements 
on a cubic lattice with 0.5-mm spacing [11]. The 
geometry was reconstructed from 1-mm sections of 
the human heart [1]. The assignment of the principal 
fiber direction was performed separately for left- 
ventricular, right-ventricular and papillary-muscle 
cells, with the fiber direction rotating from 
endocardium to epicardium [12]. Figures 2 and 3 
illustrate the amount of anatomical detail that was 
included in the model. 


Figure 1: Fiber arrangement on an yz-cross section 
approximately in the middle of the ventricles. The 
principal fiber directions are marked with arrows. 
The angles, measured from the local tangent vector 
at each element, rotate from about -50 degrees on 
the endocardium to about +45 degrees on the 
epicardium. 



Figure 2: Three-dimensional presentation of the 
ventricles as seen from the base of the heart. 

3 Results 

3.1 Simulations in slab geometry 

In order to test the propagation algorithm and 
difference approximations to calculate Eq. (1), we 
first run simulations in a three-dimensional slab (100 
x 100 x 100 cells, h = 0.05 cm). The fiber direction 
was kept constant in each simulation. Propagated 
isochrones for three cases are shown in Figure 3. 
They agree very well with theoretical predictions, 
i.e., the ratio of the longitudinal and transversal axes 
of an isochrone are proportional to the square-root 
of the conductivities in the corresponding directions. 
The propagation velocities defined from the 
isochrones yield 0.69 m/s in the axial direction, and 
0.31 m/s in the transversal direction. 


Figure 3: Propagation in slab geometry. The fiber 
direction in each case was constant, the angle from 
the x-axis being 0, 45 and 65 degrees from left to 
right, correspondingly. Each shaded zone 
corresponds to 5 ms. The ellipsoid plotted on top of 
the isochrones denotes the theoretically defined 25- 
ms isocrone. 

3.1 Conduction system 

We have created a computer interface for the 
interactive editing of the conduction system, 
utilizing the OpenGF capabilities of a SGI Visual 
PC. The organization of the conduction system was 
created by using information available from 
textbooks on human anatomy and other literature on 
the topic [13, 14]. The geometry of the conduction 
system is shown in Figure 4. The His bundle 
consists of a single branch in the right ventricle, 
whereas it resembles a fan-like sheet of fibers in the 




























left ventricle. The His bundle continues on both 
sides as the Purkinje network that contains the 
Purkinje-myocardial junction (PMJ) sites. On the 
right, the most prominent feature of the conduction 
system is the single bundle that carries the activation 
from the septum to the free wall and the papillary 
muscles along the moderator band. On the left, there 
are three major areas of activation: the septum and 
the inferior and superior free wall. There is no 
conduction network in the posterior free wall. 



Figure 4: The geometry of the conduction system 
observed from different angles. The right branch of 
the His bundle (not shown) connects the topmost 
nodes of the left and right conduction system. In the 
RV apical view (a), notice how the moderator band 
contains the conduction system activating the RV 
free wall and the papillary muscles (arrow). The RV 
free wall (b) and septal (c) views show the right 
conduction system, consisting of one main bundle 
with several smaller branches. The views of LV 
septum (d) and LV free wall (e) show the more 
complicated structure of the left conduction system. 
Notice the fan-like structure of the septal conduction 
system in (d), and the vertical spiral gap (dashed 
area) in the conduction system in (e). The PMJs are 
shown as light spheres. 

3.2 Propagation in the 3D ventricles 

The simulated activation sequence resulting caused 
by the activation of the conduction system is shown 
in Figure 5. The result agrees well with isochrones 
obtained from an isolated human heart [15]. The 
ventricular activation starts in the left ventricular 
septum (layer 110), matched by a right ventricular 
septal activation 20 ms later (layers 70-90). Almost 
simultaneously with the RV septal activation, the 
inferior (in body coordinates) and anterior LV 
activation appear (layers 90-110). These are 
followed by the activation of the RV free wall 
(layers 90-130). The RV and LV breakthroughs take 
place at 30 ms and 45 ms, respectively. In the final 
stages of the QRS, activation propagates through the 
posterior LV free wall and the pulmonary conus. 

3.3 Simulated BSPM and MCG patterns 

BSPMs were evaluated on the nodes of a 
triangulated torso model of an adult male, while the 


normal component of the cardiomagnetic field 
(MCG) was computed on a plane array in front of 
the chest. Figs. 6 and 7 show, respectively, 
simulated BSPM and MCG distributions during the 
QRS complex. In Fig. 6, the initial maximum 
resulting from left septal activation is anterior. Then 
the minimum on the back moves upward and travels 
over the right shoulder onto the right anterosuperior 
region, indicating apical activation and masking of 
the left septal activation by the corresponding right- 
septal activation. The area of positive potentials then 
drifts to the back. In the end, the activation 
propagates toward the pulmonary conus, reflecting 
as a superior positive potential in BSPM. The 
sequence of MCG patterns in Fig. 7 is also in good 
agreement with the behaviour in measured MCG 
maps. 

4 Discussion 

We created an interactive method for implementing 
the ventricular conduction system in anatomically 
accurate model of the human ventricles. The 
conduction system model was validated by 
comparing simulated isochrones, body surface 
potential maps and magnetocardiograms to existing 
data. Our validated conduction network represents 
an advancement of the forward solution of 
electrocardiography and magnetocardiography. It 
will be indispensable for simulating ventricular 
arrhythmias involving the conduction system, such 
as Purkinje-muscle re-entrant arrhythmias. The 
flexibility of our method of conduction system 
generation will allow modifications required for 
adapting our model to match measured data obtained 
during endocardial mapping. 
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Figure 5: Simulated activation isochrones. The layers are with a spacing of 10 mm, and correspond to those 
presented by Durrer et al. [15]. The isochrone spacing is 5 ms. 
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Figure 6: Simulated BSPM distributions at 10 ms steps during the QRS complex. The read and blue areas 
denote, respectively, positive and negative values. 
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Figure 7: Simulated MCG distributions at 10 ms steps during the QRS complex. The red areas denote positive 
values, where the magnetic flux enters the chest, while in blue areas the flux comes out of the chest. 
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